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ABSTRACT 

We analyzed a deep XMM-Newton observation of the radio-quiet 7 -ray PSR 
J2055-I-2539. The spectrum of the X-ray counterpart is non-thermal, with a photon 
index of r=2.36±0.14 {la confidence). We detected X-ray pulsations with a pulsed 
fraction of (25±3)% and a sinusoidal shape. Taking into account considerations on the 
7 -ray efficiency of the pulsar and on its X-ray spectrum, we can infer a pulsar distance 
ranging from 450 pc to 750 pc. We found two different nebular features associated to 
PSR J2055-I-2539 and protruding from it. The angle between the two nebular main axes 
is ~ (162.8±0.7)°. The main, brighter feature is 12’ long and <20” thick, characterized 
by an asymmetry with respect to the main axis that evolves with the distance from the 
pulsar, possibly forming a helical pattern. The secondary feature is 250” x 30”. Both 
nebulae present an almost flat brightness profile with a sudden decrease at the end. The 
nebulae can be fitted either by a power-law model or a thermal bremsstrahlung model. 

A plausible interpretation of the brighter nebula is in terms of a collimated ballistic jet. 

The secondary nebula is most likely a classical synchrotron-emitting tail. 

Subject headings: Stars: neutron — Pulsars: general — Pulsars: individual (PSR 
J2055-I-2539) — X-rays: stars 


Introduction 


The Large Area Telescope on board the Fermi satellite (Fermi-LAT lAtwood et al.ll2009l ) has 
opened a new era for pulsar astronomy by detecting 7 -ray pulsations (at E > 100 MeV) from more 
than 160 pulsars 0, ~ 25% of which are not detected at radio wavelengths. The discovery of such 
a large sample of high-energy-emitting pulsars h as enabled numerous population studies (see e.g. 
Watters &: Romanilboil I; IPierbattista et aLll 2012 L This has resulted in a better understanding of 
the physical phenomena behind the multiwavelength emission from pulsars. 
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Some of the Fermi pulsars deserved dedicated multiwavelength campaigns. Such “extreme” 
objects, accounting for the tails of the population distribution in energetics, age and magnetic field, 
are key to our understanding of the entire population, posing important limits on the different 
pulsar emission models. These campaigns have sometimes resulted in the discovery of important 
and unexpected features: e.g. the peculiar multiwavelength light curve profile of the energetic 
radio-quiet pulsar PSR J1813-1246 disagrees w ith the classical model for X-ray emission of pulsars 
( Marelli et al. 2014b : Kuiper fc Hermsen 2015 ): the isolated PSR J2021-I-4026 shows hints of long¬ 


term variations in the 7 -ray flux, challenging current high-energy emission models (jAllafort et al 

201.'ll l. 


The increase in the number of 7 -ray pulsars has been crucial to the study of phenomena related 
to highly energetic pulsars, such as X-ray and radio nebulae. Such features are usually interpreted 
within the framework of bow-shoc k, ram-pressure-dominated, pulsar wind nebulae (PWNe; for a 


review see 


Gaensler &: Slanel I 2 OO 6 I I . If the pulsar moves supersonically, shocked pulsar wind is ex¬ 


pected to flow in an elongated region downstream of the termination shock (basically, the cavity 
in the interstellar medium (ISM) created by the moving neutron star and its wind), confined by 
ram pressure. X-ray emission is due to synchrotron emission from the wind particles accelerated 


portion of the extended structure (see e.g. 

Kargaltsev & Pavlov 

2008a 1. as predicted by MHD sim- 

ulations ( 

Bucciantini 2002: Van der Swaluw 

2003: Bucciantini et al. 

20051. In order to produce a 

population of high-energy particles, a higlv-energetic pulsar is required (E greater than ~ 10^^, see 


Kargaltsev fc PavlovI l2r)n8bl ) . For classical synchrotron nebulae, we expect a relatively bright 


e.g. 

diffuse emission surrounding the pulsar, where the emissio n from the wind termination shock is 


brightest, as observed in all the o ther known cases (see e.g. iGaensler et al.l 120041 : iMcGowan et al 
20061 : iKargaltsev fc PavlovI l2008bl ) . The brightness is only slightly dependent on the interstellar 


1 /2 

medium density (oc thus we expect a relatively unifo rm brightness profile. One of the 


newly-discovered pulsars, PSR J0357-I-3205 (jAbdo et al.ll2009ll . revealed a bright tail that cannot 
be interpreted in the framework of synchrotron PWNe , due to the low energetics of the lead¬ 
ing pulsar and to a peculiar shape (jPe Luca et al.l 120111 1. Alternatively, a new model for pulsar 
nebulae was proposed, involving bremsstrahlung thermal emission from the ISM, s hocked by the 


super sonically-moving pulsar, accounting for all the characteristics of the nebula (|Marelli et al 


201,3l ). The cooldown time for this emission is orders of magnitude longer than for the synchrotron 
model, thus implying a lack of spectral variation along the tail. Pulsar bremsstrahlung nebulae 
(PEN) emission is much more dependent on the ISM density (oc n|_g^), thus pointing to a clumpy 
brightness profile, with the brightness following and enhancing the medium distribution map. For 
this type of nebula we expect that most of the energy in the pre-shock flow is carried by the 
ions, while the electron temperature is responsible for the X-ray emission: we do not expect X-ray 
emission near the pulsar, due to the Goulomb heating time of electrons. Both PWNe and PBNe 
obviously require the proper motion of the pulsar to be aligned with the main axis nebular emission, 
in the case of a co metary shape . Rece ntly, a peculiar nebular feature was found around the pulsar 
IGR J11014-6103 (iPavan et al.M201ll ). A deep Chandra observation revealed a helical structure 
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protruding from the point source, possibly associated with a jet-like emission. A smaller trail was 
detected aligned with the direction of the putative, very large (> 10 ^ km s“^) proper motion, and 
it was associated to a classical p ulsar wind nebular emission. This is reminiscent of the jet of the 
Guitar Nebula (PSR B2224-I-65, iHui &: Becked 120071 1. This feature is usually expl ained in terms o f 
a collimated ballistic jet, even if more complex explanations have been developed ([Bandieral 120081 ) . 
The helical shape and the collimated structure are the distinctive trait of this type of pulsar jet 
nebula (PJN). Another fundamental characteristic is the misalignment with the direction of proper 
motion. 


Among the extremes of the population of Fermi pulsars, PSR J2055-I-2539 (J2055 hereafter, 


Saz Parkinson et al.ll201Cll i is one of the best candidat es. R is one of the 100 brightest 7 -ray sources, 
with a 7 -ray flux of (5.5±0.1) x 10“^ erg cm“^ s“^ ( Acero et al.ll2015l ). It is located far from the 
Galactic plane, at a latitude ~-13°, making the search for counterparts at other wavelengths easier. 
A blind search using Fermi data led to the una mbiguous detectio n of the timing signature of a pulsar 
with P ~ 320 ms and P ~ 4.1 x 10“^^ s s“^ (jAbdo et al.ll2013l ). J2055 is among the less energetic 
and oldest pulsars in the Fermi zoo, with a spin-down energy of E = 5.0 x 10^^ erg s~^ and a 
characteristic age Tc = 1-2 Myr. The 7 -ray “pseudo-distance” relation (jSaz Parkinson et al.ll2010l b 
based on an empirical relation among 7 -ray pulsar luminosities an d their spin-down energies, yields 
an estimated distance of 600 pc. The second Fermi pulsar catalog (lAbdo et al.ll2013lj reports J2055 
to be one of the very few 7 -ray pulsars with a strong detection of off-pulse emission, characterized 
by a spectral cut-off: this emission probably originates from its magnetosphere, opening a new 
conundrum for magnetospheric emission models. In fact, for young pulsars this is a serious challenge 
to outer-magnetosphere models radiating only above the null-charge surface; such weakly pulsed 
emission should be rare, being expected only for nearly-aligned pulsars. 


Due to its peculiar characteristics, after two short XMM-Newton exploratory observations in 
2013 we obtained a long, uninterrupted XMM-Newton observation to study the X-ray counterpart 
of J2055 and to search for X-ray pulsations. Section [2] reports the characteristics of the observations 
and the analysis methods. In Section [3] we report the analysis of the serendipitous sources in the 
field. Section [4] describes the spectral and timing characteristics of the pulsar in the X-ray band. 
This campaign revealed the presence of spectacular nebular features, described in Section [5l Pulsar 
and nebular modeling are treated in sections [ 6 ] and [3 In the appendices, we describe the details of 
the field analysis (Appendix [A]) and of the new script we developed for our XMM-Newton analysis 
(Appendix |B|). 


2. Observations and Data Reduction 


Our deep XMM-Newton observatio n of J2055, lasting 136.2 ks, was performed on 2013 May 1 
(obs. id 0724090101). The PN camera ( Struder et al. 2001) of the EPIC instrument was operating 
in Large Window mode, with a time resolution of 47.7 ms on a 27’xl3’ field of view (FOV). The 
high time resolution combined with the large FOV allow for both the timing analysis of the J2055 
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pulsar and the spatial analysi s of the nebular structures. The Metal Oxide Semi-conductor (MOS) 
detectors (jTurner et al.ll200lh were set in Full Frame mode (2.6 s time resolution on a 15’ radius 
FOV). The thin optical filter was used both for PN and MOSs. We also analyzed XMM-Newton 
observations 0605470401 and 0605470901, taken on 2009 October 26 and on 2010 April 21, and 
lasting 24.5 and 17.9 ks respectively. The PN camera was operating in full frame mode, with a 
time resolution (73.4 ms). A timing analysis for this pulsar would be possible even using full-frame 
PN data, but it would result in a lower resolution of the light curve. We therefore used data from 
all the observations for the following analysis, with the exception of the timing analysis for which 
we considered only PN data of the longest observation. 

We used the XMM-Newton Science Analysis Software (SAS) vl3.5. For the standard analysis, 
a tool included in the SAS epproc usually randomizes the arrival time of each event within the 
time resolution. Here, we chose not to randomize the arrival times in order to achieve a better 
timing resolu t ion for our pulsar. A discussion on the implication of this choice is reported in 
Marelli et al.l (l2ni4bl ). We performed a standard analysis of high particle background with the 
SAS tool bkgoptrate (also used for the compilation of the 3XMM source catalog!^. This tool 
searches for the point at which the maximum signal-to-noise (S/N) ratio is achieved for the given 
background time series: time bins above this threshold are exclude d. We then cross-checked the 
results using an independent method, following iDe Luca et al.l (j2005l l. Both the analyses revealed a 
high contamination from flares during almost the entire observation 0605470901, that was therefore 
excluded from subsequent analysis. Moderate contamination was detected in the other observations, 
lowering the good time interval to 132 ks. Following the prescription from the 3XMM source catalog, 
for the PN camera we selected 0-4 pattern events in the 0.5-10 keV energy range and 0 pattern 
events in the 0.3-0.5 keV energy range. We applied standard flags to remove bright columns. We 
selected 0-12 pattern events for the MOS detectors in the 0.3-10 keV energy range, applying a similar 
reduction for bright columns. Finally, we excluded out-of-FOV events and residual contaminants. 
Figure [T] shows images of the counts satisfying our event selection criteria from all the datasets. 
For the timing analysis, we used the SAS tool barycen to barycenter the PN events at the best-fit 
XMM-Newton pulsar position (see Section [3|). For each spectrum, we generated ad-hoc response 
matrices and effective area files using the SAS tools rmfgen and arfgen. 


We checked for other X-ray observation of the source and found a Suzaku (jMitsuda et al.l 1200711 
observation (obs. id 405015010) taken on 2010 October 29 lasting 31.1 ks. The faintness of the X-ray 
counterpart, the crowding of the field (see Appendix E]), coupled with the low spatial resolution 
of Suzaku (the half-power diameter of the XIS instrument, onboard Suzaku, is ~ 2') make this 
observation useless for our purposes. It was thus excluded from our analysis. 


^http://xmmssc.irap.omp.eu/Catalogue/3XMM-DR5/3XMM_DR5.html 
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Fig. 1.— 0.3-10 keV sky maps showing only data used for the analysis in this paper: this includes 
the pattern, flag and proton flares filters described in Sections O [3l PN and MOSs cameras images 
from the three observations were summed. 







3. Field Analysis 


Figure [2] shows the 0.3-10 keV XMM-Newton FOY. Source detection using maximum likelihood 
fitting was done both simultaneously and on each of the EPIC-PN, MOSl, and MOS2 datasets, 
using the SAS tool edetect_chain. The resulting detections were then confirmed using the CIAO 
tool wavdetect, following prescriptions from the CIAO online cookbook [§. 


The best position of the pulsar X-ray counterpart is 20^55™48^.89 -|-25°39'57.8" (0.3"-|-1.5" 
la statistical plus systematic errors). From the most recent compilation of the LAT y-ray Pulsar 
Timing Models the best y-ray position of the pulsar is 20^55^48^.94 -|-25°39'59.1" (0.5" la 
statistical plus systematic error), fully compatibl e with the X-ray po sition. According to the logN- 
logS distribution of the Chandra galactic sources (INovara et al.ll2009ll at mid Galactic latitudes, we 
can estimate the probability of a chance detection of an X-ray source within the LAT y-ray timing 
error box and with X-ray flux comparable or greater than the X-ray counterpart (~ erg cm“^ 

s“^) to be less than 10“®. Thus, based on positional coincidence, we consider our identification 


secure. 


In Figure [2] we marked the brightest sources in the PN field of view. A bright nebular source 
is apparent in the field. This object, possibly associated with a galaxy cluster, will be described in 
Gastaldello et al. (in preparation) and not further described in this paper. We chose all the sources 
with an absorbed flux greater than 5 x 10“^® erg cm“^ s“^, with a ba detection and with at least 
150 XMM-Newton net, total counts. We also considered all the sources detected in the nebular 
regions. The study of the line-of-sight absorption of such sources could allow us to constrain the 
pulsar distance. Indeed, after selecting candidate Active Galactic Nuclei (AGN) in the FOV, it 
is possible to measure, from their spectra, the total Galactic column density in the direction of 
J2055. Next, using the pulsar column density, we can get an estimate of its distance with respect 
to the edge of the Galaxy. Such estimate could be refined if bright X-ray-emitting and optical- 
emitting stars (with known distances) were also present (see e.g. iMarelli et al.ll2013l . l2014al lbll. The 
sp atial resolution of thi s measurement is much higher than the one from the 21 cm Hj sky survey 
of iKalberla et al.l (|2005l L The resolution of the survey is 0.6°, larger than the entire XMM-Newton 
FOV: the presence of local molecular clou ds could affect the H t result even by orders of magnitude 
(see e.g. the case of PSR J1813-1246 field Marelli et al. 2014b). 


Based on the study of the spectra of serendipitous sources and of associated optical coun¬ 
terparts, we classified serendipitous sources as AGNs or candidate stars (see Appendix]^. Our 
exercise allowed us to identify 24 very likely AGN and 13 stars from a total of 58 X-ray sources. 
In particular, one of the brightest sources (source 1 in Table 1) was originally detected as a single 
source: only through the application of a new script were we able to determine that it consisted of 
the superimposition of multiple sources (see Appendix [B] for details). 


^http://cxc.harvard.edu/ciao/threads/all.html 

^https://confluence.slac.Stanford.edu/display/GLAMCOG/LAT+Gamma-ray+Pulsar+Timing+Models 
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The simultaneous fit of AGN spectra, with the column density linked, results in a relatively 
low value of column density: (2.07±0.08) xlO^^ cm“^ {la error) as expected in a mid-latitude 
region. This result is higher th an the value of 1.2 x 10^^ cm“^ obtained from the 21 cm Hi sky 
survey of iKalberla et al.l (j2005l ~). 
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Table 1:: Analysis of the brightest serendipitous sources in the XMM- 
Newton FOV. The table shows the best source position, the statistically 
acceptable spectral model fit(s), the fitted column density(s), the X- 
ray-to-optical flux ratio and the resulting identification. Typical XMM- 
Newton positional errors are 3”. We reported the la error on column 
density. 

a : Optical-to-X-ray flux ratio for the different fitting models. Where 
present, we used the V-band optical flux. If not detected, we used the 
B-band optical flux. In case of an R-band detection, we obtained the 
V-band flux (see Appendix [All 

b : Identification of these sources is based on optical and hardness ratio 
tests. 
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The Pulsar 


In order to optimize the pulsar extraction radius, we evaluated the brightness profile of the 
putative counterpart to maximize both the signal-to-noise ratio (S/N) and the number of counts. 
We found the best extraction radius to be 15^', containing a total of 1020 and 636 counts in the 
0.3-10 keV energy range in the PN and the two MOSs detectors, respectively, with a background 
contribution of ~25%. We generated ad-hoc response matrix and effective area files using the SAS 
tools rmf gen and arf gen. Due to the relatively low statistics and to optimize the fit, for the spectral 
analysis we added all the data from the MOS cameras using the HEAsoft tools mathpha, addarf 
and addrmf to sum spectra, effective area files and response matrix, respectively. PN spectra were 
analyzed separately due to the different operation modes (Large Window and Full Frame) requiring 
different response matrices. The background region was extracted from a nearby source-free region 
common for every camera and on the same CCD of the source. For the pulsar emission, we tried 
a simple power-law (pow), a simple blackbody (bb), a simple magnetized neutron star atmosphere 
(nsa), the combination of a power-law and a blackbody (bb) as well as the combination of a power- 
law and a magnetized neutron star atmosphere model (nsa in XSPEC - assuming a neutron star 
with a radius of 13 km, mass of 1.4 Mq and a surface magnetic field of 10^^ G). The simple power- 
law description gives an acceptable fit (reduced chi square x^g^=1.05, 61 degrees of freedom, dof) 
with a column density of Niy=2.18±0.26 x 10^^ cm“^ and a photon index of r=2.36±0.14 (Icr 
confidence errors). The best-fit unabsorbed 0.3-10 keV flux is (3.434 =0.27) x 10~^^ erg cm“^ s“^. 
This flux is in agreement with the one reported in lAbdo et ahl ( 20131 1 of 4 . 3 I 28 ^ 10“^^ erg cm“^ 


s“^. Figure [3] shows the spectra with the best-fit power-law model. Neither the single component 
blackbody nor the nsa model give an acceptable fit (x^g^=2.05 respectively, 61 dof) and 

were therefore rejected. Although the combined power-law plus thermal models give acceptable hts 
(both the bb and nsa give y^^^=0.99, 59 do f), the thermal normalization is compatible with 


zero 


at 2a. Moreover, an F-test (|Bevingtonlll969l ) performed comparing the simple power-law with the 


composite spectra does not point to a significant improvement by adding a thermal component. 
We thus consider the single component non-thermal model to be the best for our source. 

To search for X-ray pulsations, we used the most recent Fermi ephemeris [§. As reported in 
Section [2l we used only barycentered PN events from the 2013 XMM-Newton observation. We de- 
cided to ap ply a photon weighting technique similar to the one used for Fermi-LAT (using gtsrcprob, 


IPP 


better b ackground rejection and improves the sensitivity to pulsations, a s discussed in 
(|2ni4bl ). To this end, we further developed the Python tool reported in iMarelli et al 


discussed in 

Vlarelli et al. 

Marelli et al. 

(2014b) ("see 


Appendix IBI) . Considering only photons with a probability greater than 0.01 of comi ng from the 


pulsa r, we used a weighted Markov Chain Monte Carlo algorithm (MCMC, see e.g. I Wang et al 


2 OI 3 I ) to test and refine the 7 -ray ephemeris. Due to the shortness of the X-ray observation, we 
fixed the derivatives of frequency to ones of the 7 -ray ephemeris. We obtain a best-fitted frequency 


■"http: //f ermi. gsf c .nasa.gov/ssc/data/access/lat/ephems/ 
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Fig. 2.— Exposure-corrected sky maps showing events from all the considered datasets (left, see 
Section [2]) and only from the MOS2 datasets of the deep XMM-Newton observation (right). While 
the hrst image shows the entire dataset we used, the larger FOV of the MOS2 camera allows for a 
better view of the nebulae. Events have been extracted as described in Section [2l the images are 
in logarithmic scale. Left panel. We reported the extraction regions used for serendipitous sources 
analysis, with a color code showing the most probable association. White: the pulsar; cyan: AGN 
association ; orange: star association; green: no association. Right panel. We reported the best-fit 
main axes of the nebular features together with the best-fit length (yellow lines). We also indicate 
the areas used for the spatial analysis of the nebulae: "AMA" in green and "OMA" in yellow (see 
Section [5]) . 
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F = 3.129286±0.000003 s“^, in agreement with the 7 -ray one. With six harmonics, the result¬ 
ing H-value of 59.09 gives a tail probability of 1.25 x 10“^^. Even taking into account the ~ 50 
trials, we found the X-ray pulsation with a confidence greater than 6 cr, thus further confirming 
the identification of our X-ray candidate with the 7 -ray pulsar. Unfortunately, the low statistics 
hamper a detailed en ergy-dependent timing analysi as was done in the case of the CTA-1 pulsar 
(jCaraveo et al.l 120101 ). This could have allowed us to unambiguously accept or exclude a possible 
low-energy (e.g. thermal) component of the spectral model. 


5. Nebular emission 


The most significant feature in the XMM-Newton FOV (Figure [2|) is a thin, straight, extended 
emission passing through half of the detectors north to south and protruding from the position of 
J2055. This feature is detected in each dataset, so that we exclude an instrumental origin. Much 
less obvious is another elongated extended emission protruding from J2055 is visible south-east to 
north-west. Unfortunately, this fainter seconda ry feature is heavil y contaminated by bright point 
sources, hampering a classical analysis (see e.g. iMarelli et al.l 120131 ). 


In order to better disentangle the nebular emission from the background and the numerous 
superimposed point-like sources, we developed a new python-based tool to simulate and subtract 
their contribution from XMM-Newton data and images, based on their positions and spectra (see 
Appendix [B]). After running our tool, for each observation we created a new dataset that excludes 
photons based on the probability, P, that they come from one of the considered sources (pulsar, 
serendipitous sources and background): using a Poisson distribution, each photon has a probability 
P of being excluded. The remaining photons then come from either statistical fluctuations, or other 
sources not considered, such as our extended features. We then based our analysis of extended 
sources on this reduced dataset, subsequently confirming our results on the original datasets. 


5.1. Spatial analysis 

For each camera, we produced brightness profiles of each nebula along the main axis (AMA) 
of the feature and orthogonal to the main axis (OMA). At this point, we chose the main axis only 
through a visual evaluation, ending at the pulsar position. In order to maximize the statistical 
significance of each bin, we carefully evaluated the bin size and the width of extraction boxes, 
choosing a scale of 100” x 50” AMA and 1000” x 10” OMA for the main feature (showin in orange 
and magenta, respectively, in Figure [5]) and 100” x50” AMA and 400” x20” OMA for the secondary 
feature (shown in yellow and white, respectively, in Figure Ej). We extracted the number of counts 
in each box, summing the results of each dataset. In order to take into account the spatial evolution 
of the exposures, as well as the lack of data from an entire camera at some points, we also summed 
the exposure maps, produced through the SAS tool eexpmap, and we built a graph of the exposure 
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Fig. 3.— XMM-Newton spectra of PSR J2055+2539. The pulsar PN spectrum from the first 
observation is in black, the PN spectrum from the second observation is in red and the merged 
MOSs spectrum is in green. The best-fit spectral model is shown. 
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Fig. 4.— XMM-Newton folded, weighted light curve of PSR J2055+2539. We normalized curves 
to have the mean of the distribution equal to 1. The top curve is produced using the entire energy 
range, the middle one using only photons with energies below 1.5 keV and the bottom one using 
photons with energies greater than 1.5 keV. 









- 15 - 



Fig. 5.— Here, we show the sum of M0S2 XMM-Newton images. We applied a Gaussian smoothing 
with a kernel radius of 5”. Events are selected as described in Section [2l Upper Panel: This panel 
shows the main nebular feature as seen by XMM-Newton. Middle and Lower panels: These panels 
show the main nebular feature with the point-like sources and background removed (see Section 
[5|). The lower panel also shows the regions used for the spatial analysis in Subsection 15.11 
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Fig. 6.— Here, we show the sum of XMM-Newton images. We applied a Gaussian smoothing with 
a kernel radius of 5”. Events are selected as described in Section [2j Upper Panel: This panel 
shows the secondary nebular feature as seen by XMM-Newton. Middle and Lower panels: These 
panels show the secondary nebular feature with the point-like sources and background removed 
(see Section [5]) . The lower panel also shows the region used for the spatial analysis in Subsection 

O 
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evolution AMA and OMA, using the same boxes as before. This graph was then normalized at 
the mean value of the exposure. Finally, we divided the counts profile by the exposure profile to 
obtain a normalized brightness profile, in normalized counts arcsec“^. Then, for each nebula we 
preceded iteratively by slightly changing the main axis in order to find the maximum of the OMA 
profile compatible with zero. Therefore, we will consider this as the main axis of each nebula and 
its error as the propagated error from the maximum of the OMA profile. The brightness profile for 
both the main feature are shown in Figure [3 and [8] and for secondary feature are shown in Figure 
P and [TUI 

First, we tested the presence of the fainter nebula. To do this, we chose visually an ad-hoc 
area containing most of the nebular photons, we considered different nearby regions with the same 
shape and on the same CCDs as the secondary nebula extraction area. After considering the 
different exposure as well as instrumental effects, we obtain an expected number of counts from 
background in the secondary nebula extraction region of ~ 12600. We estimated, conservatively, 
that ~10% of the total counts are due to residual photons from point-like sources superimposed 
on the nebula, thus obtaining a total contribution of ~ 800 counts. We detected a total of 15130 
counts from the secondary nebular region (after taking into account all the instrumental effects). 
The tail probability of obtaining 15130 counts from a Poisson distribution peaked at 13400 confirms 
the presence of this feature at more than 7a. A similar calculation performed on the raw data, 
with 95% of the counts from point sources removed using circular regions, gives as well a nebular 
detection at more than 5cr. We considered the possibility that such an excess comes from point¬ 
like sources not detected by the source detection algorithms described in Section O These cannot 
realistically account for such a high number of additional counts: by taking into account the upper 
limit on the detectable flux, we expect no more than ~100 total counts from each of them, thus 
requiring a large number of unresolved sources. 

We obtain the main axis of the main nebula to be at (-0.1±0.2)° with respect to the north- 
south axis, counterclockwise. In the same reference, the secondary nebula axis is at (162.7±0.7)° 
(lu errors). 

As shown in Figure [TJ along the main axis the main nebula is characterized by a slow increase 
in flux until ~9’, followed by a rapid drop until it becomes undetectable at ~12’. Orthogonal to the 
main axis (Figure [8|), the nebular profile is clearly peaked, following a Lorentzian-like distribution 
that fades symmetrically to the axis, ~20” from it. We note that the OMA distribution is, at 
the limit, compatible with the theoretical XMM-Newton PSF distribution lldof): this 

points to a thin (no more than ~20”-thick), straight nebula or a thinner, curved nebula. To detect 
a possible main feature curvature with respect to the main axis, we extracted counts from boxes 
of 170” X125” AMA. Each box was than divided into 15” x 125” sub-boxes, following the same 
procedure as before for the exposure, then we fitted each box with a Gaussian. All the fits are 
statistically acceptable and the values (and errors) of the maximum position are reported in Figure 
[3 The resulting curve is not compatible with a constant model (Xred~^-^> 6dof): the main nebula 
shows a curvature in the clockwise direction at ~3’, then another in the counter-clockwise direction 
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at ~9’. 

Along the main axis, the secondary nebula is characterized by an almost flat flux until ~150”, 
then slowly increasing until ~250”, where it suddenly fades. Orthogonal to the main axis, it is 
peaked following a Gaussian-like distribution with a sigma of 27”±5” {la error). It does not follow 
the theoretical XMM-Newton PSF distribution (x^g^=17.3, 6dof): it is thicker (at least 50”-thick) 
than the main nebula. The low statistics prevent us from carrying out a more detailed spatial 
analysis. 

Brightness profiles of the tails in the soft (0.3-1.5 keV) and hard (1.5-6 keV) energy bands 
show no statistically significant differences from the total profiles. 


5.2. Spectral analysis 

For the spectral analysis, we evaluated the best nebular extraction regions from the results 
of the spatial analysis. For each nebula, we chose extraction boxes that ensure that the nebula is 
seen with more than a 3a significance. This results in a box of 700” x50” for the main feature and 
a box of 250” X 80” for the secondary feature, oriented as their main axes. These boxes were also 
divided in two equal parts along the main axis in order to study the nebular spectral evolution. For 
spectral analysis, we only took into account cameras that cover more than 50% of the considered 
region. The normalizations of nebular spectra were left free to vary for each camera, due to the 
different FOV and coverage of the instruments. We extracted background counts from regions with 
the same shape and angle, passing through the same CCDs as the source regions. We chose these 
ad-hoc regions to be as near as possible to the source regions, but excluding the nebular features 
and the extended feature west of the main nebula. In order to be conservative, we removed from 
the region 5”-radius circles centered at the position of point-like sources, thus excluding possible 
residuals. We followed the same procedure as in Section H] to produce spectra, but we chose not to 
add them due to the higher statistics. Spectra were grouped in order to have at least 150 counts 
per bin: owing to the high background, this allows us to have at least 25 source counts per bin in 
each camera. 

A power-law model fits well both nebulae. The best-fit absorption columns of the nebulae 
and pulsar are compatible. Then, under the hypothesis of association between the nebulae and 
the pulsar (see Section [6] for more details), we fitted all the spectra together, linking the column 
density. In this way, we find the column density of the system to be (2.22±0.17) x 10^^ cm“^ 
(90% confidence error). We note that a thermal bremsstrahlung model fits well both the nebulae, 
with column densities again in agreement with the one of the pulsar. We found no sign of spectral 
evolution with distance from the pulsar for both the nebulae down to a variation of 0.1 and 0.3 
(90% confidence) in the photon index of the main and secondary nebula, respectively. See Table 2 
for details on the fits. 


- 19 - 


AMA main nebula 



0 200 400 600 800 

Distance (arcsecs) 


Fig. 7.— Results from the spatial analysis along the main axis of the main nebula. Here, and in the 
following hgures, we plot the exposure-corrected counts versus the distance from the pulsar. Upper 
Panel: Using the original XMM-Newton datasets. Middle Panel: Subtracting the point-like sources 
and the background (see subsection 1^. Lower Panel: Difference between the best-fit maximum 
in each box and the total maximum, using a Lorentzian plus constant fit (see subsection 15.1|) . 
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Fig. 9.— Results from the spatial analysis along the main axis of the secondary nebula. Upper 
Panel: Using the original XMM-Newton datasets. Lower Panel: Subtracting the point-like sources 
and the background (see subsection [5 :td. 
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OMA secondary nebula 



Fig. 10.— Results from the spatial analysis orthogonal to the main axis of the secondary nebula; 
here, we subtracted the point-like sources and the background. We reported the best-fit Gaussian 
model. 
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src 

Nh 

r/kT 

flux“ 



10^1 cm- = 

-/keV 

erg cm“^ s“l 


total 

2.22±0.17 

- 

- 

- 

psr 

2.18±0.26 

2.36±0.14 

3.43±0.27 

1.05/61 

pwnl 

2.25±0.23 

1.82±0.08 

28.7±1.9 

1.19/249 

pwnl-br 

1.68±0.16 

^ 41+0.92 

23.5±1.8 

1.15/249 

pwn2 

9 I1-t0.73 
^•-^■^-0.59 

1.62±0.21 

5.45±1.22 

1.14/139 

pwn2-br 

J--DU_o 40 

10.3±4.4 

4.85±0.84 

1.15/139 

pwnl a 

- 

1.75±0.12 

_b 

1.11/786 

pwnlb 

- 

1.82±0.11 

- 

- 

pwn2a 

- 

1.50±0.33 

- 

- 

pwn2b 

- 

1.77±0.20 

- 

- 


Table 2: Results from spectral analysis of pulsar and nebulae. All the errors are at la confidence 
level. Here, we report the column density obtained by fitting the pulsar and the two nebular 
spectra simultaneously (total). We report the pulsar spectrum. We report the spectra of each of 
the nebulae (the main, 1, and the secondary, 2). Lastly, we report the contemporaneous spectra 
of the two parts (nearer to pulsar, a, and furthest from to pulsar, b) of the two nebulae, with the 
column density linked. 


“Unabsorbed 0.3-10 keV flux. Here, we report the flux of the instrument with the more complete view of the nebula 
(respectively MOS 2 of the second observation and PN of the second observation for the main and secondary nebula). 
*’The flux of the two parts of the nebulae is useless in order to describe them. For a complete spatial analysis see 
subsection o 
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6. Discussion 


6.1. The pulsar 

In this paper we studied the X-ray counter part of PSR J2055 -I-2539. While the majority of X- 
ray emitting isolated neutron stars (NS, see e.g. lKasoi et alJl2006l ) exhibit thermal emission coming 
from the whole NS surface, no such component is required to fit our XMM-Newton spectra. Using 
the blackbody model, for a NS radius of 10 km and a distance of 600 pc, the 3a upper limit on the 
surface temperature is 7.5 x 10^ K. This low temperature fits well the behavior of ol d pulsars (e.g. 
the J0 357-I-3205 pulsar - Morla hereafter - upper limit temperature is half this value. iMarelli et al 


([ 2013 )): we note that the cooling mechanism is highly dependent on the ma gnetic field of th e 
pulsar, so that different magneto-thermal evolution paths are expected (see e.g. iPons et al.ll2009l i. 
Taking into account the thermal radiation from hot spot(s), straight estimates of neutron star 
polar cap size based on a simple dipole magnetic field geometry, give a polar cap radius Rpc = 
R(Rn /c)^'^^, where R is th e neutron star radius, 12 is the angular frequency, and c is the speed of 
light (jPe Luca et al.ll2005l ). In the case of J2055 we obtain a value of 400 m. In that case, the 3a 
upper limit on the surface temperature is 3.0 x 10® K, too high to rule out the existence of such a 
type of emission for J2055. We note that the low-energetic J2055 is close to the death line for the 
production of (responsible for polar heating) by curvature radiation photons, thus making this 
type of emiss ion unlikely. A pul sar with comparable energetics, J0357-I-3205, shows reduced polar 
cap heating (IMarelli et al.ll2013l ). The non-thermal luminosity of J2055 is in agreement with the 
empirical re lation between the X-ray luminosity of rotation-p owered pulsars and their spin-down 
luminosity ([Possenti et al.l l2002l : iKargaltsev fc: PavlovI l2008bl i . J2055 also fits well the empirical 
relation between the X-ray luminosity (L^,) and 7 -ray luminosity (L.y) of radio-quiet p ulsars, with a 
L.t,/L .^ =1570±110 [la confidence): the mean value for radio-quiet pulsars is ~2500 (|Marelli et al 

2o3). 


We detected X-ray pulsations, thus further confirming the association of the X-ray source, with 
more than Gu significance. The shape of the pulse is consistent with a sinusoid (x^g^=1.17, 14 dof), 
although the statistics are too low for a complete characterization. The pulsed fraction is (25±3)% 
{la significance). The pulsed fractions in 0.3-1.5 keV and 1.5-10 keV energy bands are consistent, 
(28±4)% and (21±5)% respectively. A double-component pulsed emission us ually show differen t 
shapes and peak position in the thermal and non-thermal components (see e.g. iKasoi et al.ll2006l i. 
thus the consistency among pulsed fractions calls for a single-component pulsed emission, making 
the pulsation consistent with a non-therma l origin. We note that onl y ~25 pulsars in the literature 
show non-thermal X-ray pulsation (see e.g. iKuiper fc Hermsenll2ni5l ). fewer than half of these have 
a 7 -ray counterpart and only three of them are radio-quiet (Geminga, Morla and J1813-1246). This 
makes J2055 an interesting object to test models for high-energy emissions of pulsars. 
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6.2. The nebulae 


Most of the analysis we carried out focused on the two elongated ragions of extended emission 
we detected inside the XMM-Newton FOV. Both the sources have morphologies typical of PWNe 
or jets from pulsars, with an elongated shape, brighter at one edge and almost symmetrical with 
respect to the main axis. Their spectra are non-thermal, best fitted by a power-law, thus confirming 
their association either with nebulae from pulsars or a jet from an AGN. An interpretation of the 
feature as an AGN jet can be safely discarded, owing to the angular extent of the feature: it 
would imply an unrealistic physical size, unless the source is quite local (a huge 200 kpc long jet 
would imply an angular scale distance of order 40 Mpc for the brighter feature and 150 for the 
fainter, assuming standard cosmological parameters), which would call for a rich multiwavelength 
phenomenology (the host galaxy itself - with an angular scale well in excess of 1 arcmin - should be 
clearly resolved by ground-based optical images we used). We thus refer to these features as pulsar 
nebulae. 


The morphology of the nebulae, apparently protruding from J2055 and connected to the pulsar 
counterpart, strongly argues for a physical association of the systems. Three other sources, la, lb 
and 2 in Table 1, can call for such a morphological association for the fainter nebula. Thanks to 
accurate spectral and multiwavelength analysis, we concluded that source lb and 2 are stars so 
that they cannot be associated to the features. The spectrum of source la, on the other hand, 
is non-thermal and we found no optical counterpart that can be associated, up to a magnitude 
of 21. These characteristics point either to a pulsar or an AGN. Taking into account the low 
number of known X-ray pulsars (to date, fewer than 200) in the entire XMM-Newton sky, the 
probability of finding a serendipitous pulsar within 1’ of J2055 is negligible, thus confirming the 
association of both nebulae with J2055. Moreover, the best fitted column density of source la is 
higher than the Galactic column density we found ((2.07±0.13) xlO^^ cm~^), thus pointing to an 
extragalatic object. As further confirmation of the nebulae-pulsar association, their column density 
are in agreement, 2.18±0.26 x 10^^ cm~^, 2.25±0.23 x 10^^ cm~^ and 2 . 11 ^Q 5 g x 10^^ cm“^ (90% 
confidence) respectively for the pulsar, main nebula and secondary nebula. 


_ The Galactic absorption column in the direction of the pulsar, predicted bv iDickev X Lockman 

( 19901 ) at ~450 pc and based on the Hj distribution (1.2 x 10^^ cm~^), is lower than the best-fit value 
for the pulsar-nebulae system ((2.22±0.17) xlO^^ cm~^), pointing to a lower limit for the pulsar 
distance of ~450 p c. This result is in broad ag reement with the value estimated by scaling the 7 -ray 
flux of the pulsar (ISaz Parkinson et al.ll201fll ): the method hinges upon the observed correlation 


betw e en the intrinsic 7 -ray luminosity and spin-down power of the pulsar (jSaz Parkinson et al 


2010 ; 


Abdoetah 


2 OI 3 I). assuming a beam correction factor of one for the 7 -ray emission cone of 


all the pulsars (IWatters et al.ll2009l i. By applying this relation, we have a 7 -ray efficiency of 0.04, 
an X-ray efficiency of 3 x 10““^ and a distance of 600 pc. We obtain an upper limit of ~750 pc by 
requiring the 7 -ray efficiency to be less than 1 . Another important piece of information could come 
from the serendipitous stars we found in the field: their distance estimates from optical analyses 
could provide limits on the distance of the pulsar. Future multi-filter optical observations will 
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provide this limit. 


We found two nebulae associated with J2055, with angular dimensions of 12’ x 20” and 250” 
X 30”, respectively, for the brightest and faintest. The observed extensions of the nebulae, at a 
distance of 600 pc, would correspond to physical dimensions of ~ (2.1 x 0.05) pc and ~ (0.7 x 0.09) 
pc respectively (assuming no inclination with respect to the plane of the sky). The luminosities 
of the nebulae in the 0.3-10 keV energy range (assuming d = 600 pc) are 1.2 x 10^^ erg s“^ and 
2.4 X IqSO 

erg s corresponding to fractions of 2.4 x 10 ^ and 3.0 x 10 ^ of the pulsar spin- 
down luminosity. A few elongated tails of X-ray emiss i on associated to rotat i on-powered pulsars 
have been discovered in the past (jGaensler et al.ll2004l : iMcGowan et al.ll2006l : iBecker et al.ll2006l : 
Kargaltsev &: Pavlov! l2008bl b Although for our pulsar we have no information about the proper 
motion, the bow-shock PWN scenario would seem the most natural explanation for one of the 
nebulae. Luminosity values are fully compatible with that measured for other synchrotron nebulae, 
for which the pulsar channel into their tails 10“^ to 10““^ of their rotational energy loss. In the 
case of synchrotron emission tails, if we assume the optimistic maximum Lorentz factor of ~ 10® 
for electron acceleration in such a low-energetic pulsar magnet osphere, and the hig h value of the 
ambient magnetic field of ~ 50 /rG (following considerations in iDe Luca et al.ll2011h . it is possible 
to estimate the synchrotron cooling time of the emitting electrons as Tgync ~ 100 (B/50/iG)“®/^ 
(E/lkeV)“^/^ yr. Coupling this value with the estimated physical length of the feature yields an 
estimate of the bulk flow speed of the emitting particles of 20,000 and 3,000 km s assuming 
no inclination with respect to the plane of the sky. The first value is only marginally consisten t 
with results in literature, while the second one is fully consistent (iKargaltsev fc Pavlov! l2008bl ) . 
Taking into account the energetics, both the nebulae are at least marginally consistent with a 
classical synchrotron nebula explanation. For classical synchrotron nebulae, we expect a relatively 
bright diffuse emission surro unding the pulsar , where the emission from the wind termination shock 
is brightest. The model of iGaensler fc Sland (120061 ) predicts a low-scale termination shock of ~ 
0.6” fo J2055 (at 600 pc), thus not resolved by XMM-Newton, even in the case of typical ambient 
densities (0.01 atoms cm“®) and pulsar velocities (some hundreds km s“^). In that case, part of 
the flu x we assign to the pulsar comes instead from the termination shock (jKargaltsev &: Pavlov 
2008al f . Along their main axis, the brightness profile of both the nebulae are consistent with an 
almost flat behavior, with a sudden decrease at the end. This is acceptable for synchrotron emission 
nebulae: in fact, the pulsar wind is more energetic in the surrounding of the pulsar, where the loss 
of energy through synchrotron emission is higher, decreasing in flux and energy along the pulsar 
trail. Taking also into account an inhomogeneous medium and/or r nagnetic field, this behavior is 
expected, and observed for many objects in the literature (see e.g. IKargaltsev fc Pavlov! 120081)1 ) . 
Such a synchrotron cooling of the particles injected at the termination shock induces a significant 
softening of the emission spectrum as a function of the distance from the pulsar in bow-shock PWNe. 
Taking into account the upper limit variation of the photon index we found, for the brightest tail 
of J2055 we do not have the predicted spectral variation. Glassical theories of ram pressure particle 
confinement in pulsar tails are hardly put to the test by parsec-long nebulae, that would require 
much higher efficiencies than in all the other cases. The tightness of the main feature of J2055, if 
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confirmed to be a synchrotron nebula from its shape and the pulsar motion, cannot be accounted 
for by any model in the literature. We conclude that the low energetics of the pulsar, the lack 
of any spatial-spectral variability and the tightness of the brighter nebula disfavor the classical 
synchrotron emission nebula. This model can explain the characteristics of the secondary feature. 

If we consider that the main nebular emission comes from thermal bremsstrahlung, we can 
assume an op timistic cone angle of the tail of 1°. Taking into account the spectral best fit values 
and following iMarelli et al.l (|2013l i. this implies an interstellar medium temperature of ^ 50,000 K, 
depending on the inclination of the nebula. The pulsar velocity along the plane of the sky would 
be ~ 2300 km s“^, one of the highest values for a pulsar. The high temperature of the pre-shock 
i s consistent with that of the hot phase of the ISM, which fills a large fraction of the Galaxy 
(jBland-Hawthorn fc EevnoldsI 120001 ). The variation of the symmetry with respect to the main axis 
of the nebula with the distance from the pulsar can be explained in terms of a clumpy distribution 
of the interstellar medium. The long expected cooling time (~ 10^ yr) also explains the lack of 
spatial-spectral variation of the nebula. In any case, we need the lack of emission surrounding 
the pulsar for the electrons must be heated by the ions, heated by the pre-shock flow. We have 
no indication about the interstellar medium density, that would depend on the inclination of the 
system, but taking into account values similar to the ones of PSR J0357-I-32, we expect that within 
^ 50” of the pulsar there should be no nebular emission. The brightness profile along the main 
axis clearly rules out this model both for the main and the secondary nebulae. 

Both scenarios we considered need the proper motion of the pulsar to be aligned with the 
main axis of the nebula. The presence of two different tail-like features implies that at least one of 
them cannot be explained in terms of classical nebular models. A magnificent example of a long, 
collimated X-ray nebular feature, misaligned wi th the pulsar proper motion is the Guitar Nebula 
(the radio-loud PSR B2224-I-65, located at 1 kpc, Hui &: Becker (2007)). The Guitar feature is much 
shorter (~2’ long) and larger than the J2055 one, with the similar nebula- pulsar flux ratio of ~6. 
The main feature of J2055 is thus compatible with a Guitar-like feature from a nearer object (~150 
pc rescaling Guitar feature). This feature is usually expla ined in terms o f a collimated ballistic jet, 
even if more complex explanations have been developed (jBandieral l2008l ) . These models predict a 
spectral change across the fea ture, excluded from XMM-Newton observations in the case of J2055. 
Recently, Pavan et al. ( 2014al ft]) found a long helical jet-like structure protruding from the energetic 
pulsar IGR J110I4-6103. Two different features protrude from the point-like object, one shorter 
(1.2’) and wider probably aligned with the proper motion, also radio-emitting, and one 5.5’-long 
forming an angle of ~ 104° with the shorter feature axis. The jet scenario is the most probable, 
mainly due to the helicoidal shape of the nebula, revealed by the Ghandra observation. While 
similar to the J2055 feature in length, the fluxes of the IGR object and its two main features are 
comparable. With J2055 we are facing a much more collimated nebula protruding from a much 
less energetic pulsar, theoretically difficult to account. In order to confirm the helical pattern of 
the main nebula and proceed with theoretical modeling, as well as to confirm the misalignment of 
the proper motion, high-resolution images of the nebula should be taken at different epochs. 
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We recently obtained two new X-ray observations of the J2055 system, taken by Chandra. The 
system will be observed during a long, uninterrupted observation, allowing for a full characterization 
of the tails with a much better spatial resolution than XMM-Newton. With such an observation, 
we will also improve the results of our script on XMM-Newton data. After two years, Chandra will 
re-observe the J2055 system in a short observation. This will allow us to detect the proper motion 
down to ~ 350 km s“^, for a distance of 600 pc and a null inclination with respect to the plane of 
the sky. Such information will be crucial to understand the origin of the two nebular features of 
PSR J2055+2539. 


7. Conclusions 

We analyzed the XMM-Newton observations of the radio-quiet y-ray PSR J2055-I-2539. We 
confirmed the X-ray counterpart association of the pulsar. The X-ray emission is non-thermal, with 
an X-ray efficiency of Vx = ^xlE ~ 3 x 10“^ (for a distance of 600 pc), typical of X-ray pulsars. 
The 7 -to-X-ray luminosity ratio is ~ 1600, typical of radio-quiet pulsars. The lack of cooling 

thermal emission is due to the old age of the pulsar, while the lack of hot spot thermal emission 
can be due to its low energetics. We detected non-thermal X-ray pulsations, following a sinusoidal 
prohle with a pulsed fraction of ~ 25%. Taking into account considerations on the y-ray efficiency 
of the pulsar and on its X-ray spectrum, we can infer a pulsar distance ranging from 450 pc to 750 
pc. 

We developed a new script to simulate the spatial distribution of XMM-Newton counts from 
input point-like and extended sources Using this, we found and analyzed two different nebular 
features associated to J2055. Both nebulae have an elongated profile, protruding from the pulsar 
and almost symmetrical to a main axis. The main axis of the two nebulae are separated by an 
angle of ~ 163°. The main, brighter feature is 12’-long, <20”-thick (corresponding to 2.1 x 0.05 
pc at a 600 pc distance), characterized by an asymmetry with respect to the main axis evolving 
with the distance from the pulsar, possibly forming helical pattern. The secondary feature is 250” 
X 30” (that corresponds to 0.7 x 0.09 pc at a 600 pc distance). Both nebulae present an almost 
flat brightness prohle, with a sudden decrease at the end. The nebulae can be htted either by a 
power-law model or a thermal bremsstrahlung model. Their efficiencies (Lx/E) are respectively 
2.4 X 10“^ and 3.0 x 10“^ for the main and the secondary nebula. A plausible interpretation of 
the brighter nebula is in terms of a collimated ballistic jet, mainly due to its thickness and shape. 
The secondary nebula is most probably a classical synchrotron-emitting tail. Future Chandra 
observations will allow us to better study the shape of the nebulae and to hnd the proper motion 
of the pulsar, thus allowing us to test our interpretations. Finally, we found and studied a bright, 
serendipitous galaxy cluster with a central AGN. 
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A. Appendix A: field analysis 


Based on the study of spectra and possible optical counterparts, we can classify serendipitous 
sources as AGN or candidate stars, allowing us to constrain the pulsar distance. Indeed, after 
selecting candidate AGN in the FOV, it is possible to measure from their spectra the total Galactic 
column density in the direction of J2055. 


For this analysis, we decided to use only brighter sources, those from which it is possible to 
extract and fit a spectrum: we therefore require that they have an absorbed flux greater than 5 
X 10“^® erg cm“^ s“^, with a 5a detection and with at least 150 XMM-Newton net, total counts. 
We also considered all the sources detected in the nebular regions. As described in Section [3l we 
found 57 serendipitous bright X-ray so urces in the held of v iew. In order to classify them, we 
performed a spectral analysis, as done in iMarelli et al.l (j2014bl l. All the sources, with the exception 
of source 1, are compatible with a point-like distribution. As described in Appendix O a deeper 
analysis of source 1 revealed It as the superimposition of two or three different sources. Extended 
so urces are discu s sed e lsewhere. As a first method of classification, we used their spectra. Based 


on iNovara et al.l (120091 ) . at mid-Galactic latitude we expect that ~75% of serendipitous sources 
with the considered flux are AGN and the remaining are stars. X-ray emission from AGN is best 
described by a power-law, while X-ray emission from stellar coronae by single or double apec 
(emission spectrum from collisionally-ionized diffuse gas). For sources located inside extended 
emissions, we extrapolated the spectrum and flux of the extended source and froze it into the 
spectral fit. Out of these 58 sources, 30 can be realistically fitted only be a power-law and 5 only 
by apec(s); for the remaining 23 objects, this type of analysis is not conclusive (see Table 1). 


Since the count statistics of some of the selected X-ray sources are too low to discriminate the 
spectral model, we performed a qualitative spectral analysis using the count rate (GR). The GR is 
measured in three energy ranges: soft (0.3 — 1 keV), medium (1 — 2 keV) and hard (2 — 10 keV). 
We then computed the two different hardness ratios: 


HR12 = [CR{1 - 2) - CR(0.3 - 1)]/[CR{1 - 2) + CR(0.3 - 1)] 
R'R23 = [CR{2 - 10) - CR{1 - 2)]/[CR{2 - 10) + CR{1 - 2)] 


Large/small HR12 values points to a large/small absorption, a large/small HR23 indicates 
a hard/soft spectrum. The values of the hardness ratios are compared to a matrix of simulated 
values for apec (kT ranging from 0.5 to 5.5keV) and power-law (F ranging from 1.5 to 2.5) mod- 
els. Each spectral model is computed using the average interstellar medium absorption given by 


Dickey &: LockmanI (|l990l l. We extracted and simulated count rates from the PN camera of the sec¬ 


ond observation and the sum of MOSs cameras, then we combined them. These values are plotted 
in Eigure [11] together with a sample of the measured HR ratios. From the position of the latter with 
respect to the former, an interpolation is possible over the two models and the values of the model 
parameters. Although this method proved best for the case of the highly-absorbed serendipitous 
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sources in the field of PSR J1813-1246 ( Marelli et alJl2014bl l. for J2055 field this method did not 
improve the spectral results already obtained for most of the sources. Only in the case of source 2 
do we have an indication of a star-like spectrum and for sources 3, 43 and 59 of AGN-like. 


Another common way to confirm X-ray classification of sources is based on multiwavelength 
an alysis: the X-to-optical flux ratio is a good indicator of the nature of X-ray emitters. According 


to 


La Palombara et al.l (|2006l i , AGN have typical logarithms of X-to-optical flux ratios higher than 


-0.15(-|-0.2), while stars lower than -|-1.0(-|-0.65) in the V(B) band. Inside each X- ray source error 


box, we looked for association with optical sources from the NOMAD catalogue (jZacharias et al 


20051 1. considering the V-band magnitude as reference or B-band where necessary. Alternatively, 


we used the R-band magnitude: V magnitude was then evaluated from R using the following 
method. We considered the best-fitting X-ray models to estimate the column densities for each 
source. In case of multiple fitting models, we repeated the follow i ng fo r each model. We de - 
absorbed R magnitudes by using coefficients from iPredehl &: SchmittI (|l995l l ; iGardelli et al.l (|l989l l , 
obtaining the R-band absolute magnitude. Then, assuming a thermal spectrum and/or a power-law 
spectrum, based on the X-ray model, we used the on-line NIGMOS tool 0 to evaluate the V-band 
absolute magnitude. Few of our AGN-like objects have optical counterparts inside their error box 
that could be due to s purious matches. I n ord e r to estimate the num ber of spurious matches, we 
used the relation from ISevergnini et al.l (|2005l l ; iNovara et al.l (j2009l l . This yielded a probability 
of chance coincidence of 21%, wich means that, at our limiting magnitudes, contamination effects 
cannot be ignored. Each of the objects revealed as star-like from spectral or HR analysis has an 
optical counterpart that agrees with the expected X-ray-to-optical flux ratio. Table 1 reports the 
associated optical counterparts and expected upper limits for AGN. 


Finally, we combined results from spectral and multiwavelength analysis to classify our serendip¬ 
itous sources, with the restrictions reported above. All the sources with spectral model or HR 
compatible with a star coronal emission and with a star-like X-to-optical flux ratio are considered 
as stars: these are 13 objects. All the sources with spectral model or HR compatible with an AGN 
emission and with a AGN-like X-to-optical flux ratio are considered as stars: these are 24 objects. 
21 sources cannot be uniquely classified. 


B. Appendix B: a script for X-ray photons weighting 

We developed a script to simulate the spatial distribution of photons from input point-like and 
extended sources. After this simulation, by setting the value of the variable ‘other’, it is possible 
either i) to search for new sources using the reduced dataset without photons from considered 
sources (other=yes) or ii) to extract photons from a single source to perform high-level spectral 
and timing analysis (other=no). We note that this type of approach is particularly well suited 


c 


http://www.stsci.edu/hst/nicmos/tools/conversion_form.html 
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Fig. 11.— Distribution of HR12 vs. HR23 of a sample of selected X-ray sources. Error bars are 
reported at la. Blue stars indicate the expected HR12 vs. HR23 computed for power-law spectra 
with r from 1.5 and 2.5. Red stars indicate the expected HR12 vs. HR23 computed for apec spectra 
with kT from 0.5 to 5.5 keV. Each spectral mod el is computed using the average interstellar medium 
absorption given bv lDickev &: LockmanI (jlOOOl ). 
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to combine the high spatial resolution of Chandra-ACIS observatory with the high spectral and 
timing resolutions of XMM-Newton. The first satellite can provide the number and exact position 
of point-like sources as well as the flux and shape of nebular features: these can be used as inputs 
of the script to perform high-level XMM-Newton analysis. In this paper, first we used the script to 
clean out the nebular features from point-like sources reported in Table 1, see Figure [T2j Then, we 
used again the script to extract photons from the pulsar to perform the timing analysis. 

We used python v2.7.5 language to build the script. This python script recalls standard 
python modules - os, sys, math, shutil, multiprocessing, subprocess - as well as modules 
dedicated to scientific analysis - xspec, numpy, pyfits. This script has a number of variable 
input parameters: 

- the type of analysis: to search for new sources or to extract photons from a source (‘other’ 
variable); 

- the energy range and energy resolution to be used, in eV; 

- the spatial accuracy to be used for point-like PSF reconstruction; if one of the sources is extended, 
it should be set also the extended PSF accuracy; 

- the names and positions of input and output files. 

The script takes a list of inputs from an ASCII file: 

- the position of each source, both in ecliptic and physical coordinates; 

- the spectral model and parameters of each source, readable by xspec module; 

- the source extraction circle radius from which each spectrum of point-like source is extracted; 

- the spectrum of background, readable by the xspec module; 

- the extraction circle radius of background spectrum. 

Finally, as input it needs some FITS files: 

- event file, filtered by the user for instrumental effects and energy; 

- image file constructed from the event file; 

- exposure file constructed from the event and image files; 

- calibration file (for XMM-Newton^ the CCF file); 

- the spectrum, background spectrum, response matrix (rmf) and effective area (arf) file from each 
source and from the background; 

- if extended sources are loaded, for each of them an image of the theoretical spatial shape of the 
nebula (before the PSF application), possibly with different normalizations of each pixel. 

As output, it produces: 

- for each input source, an ASCII file containing all the commands ran by the script; 

- for each input source and the background, a column is added to the hts file: this show the 
probability (0<P<1) of the photon to come from that source; 

- for each input source, the image of the simulation of that source, see Figure [121 

- an image containing all the simulated sources, an image with residuals and an image containing 
the deviation of simulated map from counts map, in sigmas. Due to the high computing time 
needed, each source is treated separately on a different CPU and the results are merged only at 
the end, to compute probabilities and create images. 
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The first part of the script produces the mean count rate over the entire observation of each 
source, in each energy band. Here, the spectrum, background spectrum, rmf and arf are loaded 
using the pyxspec library and the input spectral model is assigned (but not fitted). Based on the 
spectral model interpolated with rmf and arf, the script computes the count rate. We note that 
these values are corrected for PSF tails and exposure map. 


The second part of the script constructs the PSF of each sour ce, in each energy band. We used 
the last theoretical model for XMM-Newton PSF, as reported in iRead et al.l ([20111) and following 
calibration documents. The normalization of the PSF function is set in order to obtain the correct 
number of expected count rates, as calculated in the first part of the script. In case of an extended 
source, this is treated as a number of point-like sources with positions and normalizations following 
the input model image. We note that, in order to reduce the computing time, a different, lower- 
resolution, spatial grid can be set for extended sources. The use of this complicated, spatial and 
energy-dependent, PSF function rather than the simple King function previously adopted resulted 
in a significant improvement in the positional accuracy, as well as spurious sources identification 
for automated XMM-Newton searches of serendipitous sources. This change in the PSF becomes 
more important for bright sources in the field, allowing for a better disentangling photons from each 
source, as well as the detection of nearby faint sources. In the case of the J2055 field, the adoption 
of this newer PSF resulted in the clear identification of the two, or possibly three, different sources 
responsible for the emission from source 1 in Table 1. 


The third part of the script computes the number of observed counts, from each source in 
each energy band, in each pixel of a spatial grid, with a step set by the user. It is based on the 
theoretical PSF function and values obtained in the second step. Finally, the obtained matrix is 
multiplied for the input exposure map to obtain the number of expected counts. 


The fourth part of the script creates output images. The selected spatial grid is collapsed to 
agree with the input counts map resolution and results from all the energy bands are summed. For 
each source, a simulated image is written and they are all summed to write the total simulated 
sky. The script also creates a residual image, subtracting the simulated map (S) from the counts 
map (C), and a significance map defined for each i bin as SIi=||Oj-Si||/\/Oi where Oj!=0 and Sli=0 
otherwise. 


The fifth part of the script computes the probabilities. The script assigns to each photon in 
the FITS file a list of predicted counts values for each source, based on the pixel of the spatial grid 
where the photon falls. If ‘other’=yes, it also takes into account the observed number of counts 
in that pixel. If higher than predicted, it adds the value (number of observed counts)-(number 
of predicted counts) to the list, otherwise it adds a zero. Finally, the sum of the values for each 
photon is normalized to 1 and the obtained list of probabilities is added to the input FITS file. 
We note that the ‘other=yes’ mode can be used only to search for new sources but not for the 
analysis of an already-detected source: in fact, this would induce a systematic error due to the 
wrong characterization of the statistics. On the other hand, accurate analysis of the residuals 
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would allow one to find low-statistic point-like sources or extended sources that can be added to 
the input list to iteratively improve the results of the script. In order to extract photons from the 
source of interest, only sources ~ 7.5’ away from it should be considered. XMM-Newton theoretical 
PSF is null further than 5’ from the source, so that only sources that affect this circle are necessary. 


In order to validate this method, we applied it to several different pulsar fields: we used the 
weighted MCMC analysis presented in Section|4]to compare the H-value resulting from weighted and 
unweighted tests. We obtained improvements of H-value ranging from 10 to 50%, with better results 
for pulsars with multi-component spectra, su perimposed on nebul ar and/or point-like sources. In 
the case of the already-published J1813-1246 (jMarelli et al.ll2014bl i. a simpler version of this script 
increased the H-value from 11123 to 12092; the use of this new method further increased it to 
12432. In the case of J2055, the H-value increased from 52 to 66 (with the same number of trials). 
In the case of Geminga, we used some of the public XMM-Newton observations to detect the non- 
thermal pulsation, in the 2-10 keV band. Using our new script, we found an increase in the H-value 
from 120 to 220, almost doubling it. Tests on other pulsars are in progress. We also considered the 
distribution of residuals of the considered point-like sources. For each of the ~ 50 object considered, 
we found residuals of no more than a few percent of the original fluxes, thus confirming the validity 
of the method. 


The case of source 1 (see Table 1) is representative of the capacity of this script. That object 
was revealed both by the SAS tool edetect_chain and the CIAO tool wavdetect as a single 
point-like emission. A brightness profile analysis does not reveal significant deviations from the 
theoretical XMM-Newton PSF. A spectral analysis gives instead a peculiar result: this source 
cannot be fitted by a single-model spectrum but requires a double-component spectrum (apec+apec 
or apec+power-law). Our script allowed for a more accurate analysis of the source. Taking into 
account a single point-like object as source 1, we obtain residuals of more than 40% of source counts 
(while the other ~50 sources give only a few percentage points), distributed in pixels south-east 
to north-west from the best source 1 position; we also obtain an overestimation of counts in pixel 
south-west to nort-east, so that the total number of simulated counts is in agreement with the 
observed ones. We note that the nearby nebula should not significantly contribute to the source 
counts, due to its faintness. Then, we tried a model with two different sources: the positions were 
iteratively refined in order to minimize the residuals. We found that two sources separated by 
describe adequately the ‘source 1 system’, with lower residuals of ~15%. The brightest one has a 
flux ~10 times larger than the weakest and its spectrum is best described by a power-law while 
the other one spectrum is best described by an apec; spectral parameters and optical counterparts 
allow us to associate them with an AGN and a star, respectively. We note that a three-sources 
system, forming an equilateral triangle with a side of ~ 10" would lower the residuals to ~3%. 


Althou gh this program is m uch more evolved than the script we used for the analysis of PSR 
J1813-1246 (iMarelli et al.ll2ni4bl ). it still needs further technical improvement to reduce the high 
computing time. A future version of this script could also be more user-friendly: a simplified 
version of the analysis could be run using the energy bands and associated count rates for each 








- 35 - 


of the objects reported in the 3XMM catalog of serendipitous XMM-Newton sources 0. Future 
implementations could also all ow the user to analy ze data from X-ray observatories different than 
XMM-Newton, e.g. NuSTAR ([Harrison et al.ll2013l h by simply changing the PSF model and some 
variables. 


This work was supported by the ASI-INAF contract 1/037/12/0, art.22 L.240/2010 for the 
project "Calibrazione ed Analisi del satellite NuSTAR”. 

Facilities: XMM. 
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Fig. 12.— Left Panel: sum of all the 0.3-10 keV XMM-Newton counts images, represented in 
logarithmic scale. Here, we did not apply all the filters described in Section [2] to remove bright 
columns; Right Panel: sum of all the 0.3-10 keV XMM-Newton images simulated by our script, 
represented in logarithmic scale. The asymmetry of the PSF is apparent for sources with high 
off-axis angles and hints of spikes are visible for the brightest sources; examples of the differences 
between the two figures are showed in Figures [5l [6] 







